﻿using System;
using System.Numerics;

namespace Algorithm_Comparison
{
    /// <summary>
    /// Runs Enhanced Stochastic Evolutionary algorithm based on the PhiP criterion on the given LHD
    /// MaxCycles is redundant, nowhere in the experiments is maxCycles used.
    /// </summary>
    class EnhancedStochasticEvolutionaryPhiP
    {
        // All the variables set to the values described in the paper
        public static readonly double SmallPercentage = 0.1;
        public static readonly double LargePercentage = 0.8;
        public static readonly double Alpha1 = 0.8;
        public static readonly double Alpha2 = 0.9;
        public static readonly double Alpha3 = 0.7;
        public static readonly int MaxJ = 50;
        public static readonly int MaxM = 100;

        public readonly LHD InitialLHD;
        public readonly int PValue;
        public readonly long MaxCycles;
        public long MaxExchanges;
        public readonly bool Euclidian;

        public readonly int J;
        public readonly int M; // #iterations in inner loop
        public readonly double InitialTreshold;

        public long CompletedCycles; // redundant
        public long NrOfExchanges;

        public LHD BestRecordedMaximinLHD;
        public LHD EndResult;
        // if maxExchanges is not 0, maxCycles is no longer used as a stopping criterium
        public EnhancedStochasticEvolutionaryPhiP(LHD initialLHD, int pValue, long maxCycles, long maxExchanges=0, bool euclidian = true)
        {
            InitialLHD = initialLHD;
            Euclidian = euclidian;
            PValue = pValue;
            MaxCycles = maxCycles;
            MaxExchanges = maxExchanges;

            InitialLHD.CalculateDistanceMatrix(Euclidian);
            InitialLHD.CalculatePointMinDistance();
            InitialLHD.CalculatePhiP(PValue);

            J = CalculateJ();
            M = CalculateM();
            InitialTreshold = CalculateInitialTreshold();

            CompletedCycles = 0;
            NrOfExchanges = 0;
            BestRecordedMaximinLHD = InitialLHD;
            Run();

        }

        /// <summary>
        /// Calculates J: the number of distinct designs to generate in each iteration
        /// </summary>
        /// <returns></returns>
        private int CalculateJ()
        {
            BigInteger tempResult = 1;

            for (int i = InitialLHD.NumberOfPoints; i > InitialLHD.NumberOfPoints - InitialLHD.Dimension; i--)
                tempResult *= i;

            for (int i = 1; i <= InitialLHD.Dimension; i++)
                tempResult /= i;

            tempResult /= 5;

            return (tempResult > MaxJ) ? MaxJ : (int)tempResult;
        }

        /// <summary>
        /// Calculates M: the number of iterations in inner loop
        /// </summary>
        /// <returns></returns>
        private int CalculateM()
        {
            return 10*InitialLHD.Dimension; // !!! Assuming Nl = Ne !!!
            //Redundant code:
            /*
            int result = (int)((2*InitialLHD.NumberOfPoints*InitialLHD.Dimension)/(double)J);

            if (result >= MaxM)
                return MaxM;

            return result;*/
        }
        
        //Calculates the initial treshold such as described in the paper
        private double CalculateInitialTreshold()
        {
            return 0.005d * InitialLHD.LastCalculatedPhiP;
        }

        private readonly static Random RandomGenerator = new Random();

        public void Run()
        {
            //Initialize
            LHD lhdX = InitialLHD;
            LHD lhdXBest = lhdX;
            double treshold = InitialTreshold;
            do
            {
                LHD lhdXOldBest = lhdXBest;
                int nAccepted = 0;
                int nImproved = 0;
                for (int i = 0; i < M; i++)
                {
                    int column = i%InitialLHD.Dimension;
                    LHD lhdXTry = null;

                    for (int jCounter = 0; jCounter < J; jCounter++) // j random element exchanges within column i mod m
                    {
                        NrOfExchanges++;
                        int[] pointIndex = new int[2];

                        pointIndex[0] = RandomGenerator.Next(0, InitialLHD.NumberOfPoints);
                        pointIndex[1] = RandomGenerator.Next(0, InitialLHD.NumberOfPoints - 1);
                        if (pointIndex[1] == pointIndex[0])
                            pointIndex[1]++;
                        LHD temp = lhdX.LocalMove(pointIndex[0], pointIndex[1], column, Euclidian);
                        temp.CalculatePhiP(PValue);

                        if (lhdXTry == null || temp.LastCalculatedPhiP < lhdXTry.LastCalculatedPhiP)
                        {
                            lhdXTry = temp;
                        }
                        if (temp.MinDistance > BestRecordedMaximinLHD.MinDistance)
                            BestRecordedMaximinLHD = temp;
                    }

                    if (lhdXTry.LastCalculatedPhiP - lhdX.LastCalculatedPhiP < treshold*RandomGenerator.NextDouble())
                    {
                        lhdX = lhdXTry;
                        nAccepted++;
                        if (lhdX.LastCalculatedPhiP < lhdXBest.LastCalculatedPhiP)
                        {
                            lhdXBest = lhdX;
                            nImproved++;
                        }
                    }
                }


                double acceptanceRatio = nAccepted / (double)M;
                double improvementRatio = nImproved / (double)M;
                if (lhdXOldBest.LastCalculatedPhiP > lhdXBest.LastCalculatedPhiP) //flagimp = 1, improvement process
                {
                    if (acceptanceRatio > SmallPercentage && improvementRatio < acceptanceRatio)
                        treshold *= Alpha1; // decrease treshold
                    else if (acceptanceRatio > SmallPercentage && acceptanceRatio == improvementRatio)
                        treshold = treshold; //redundant code, maintain treshold
                    else
                        treshold /= Alpha1; // increase treshold
                }
                else//flagimp = 0, exploration process
                {
                    if (acceptanceRatio < SmallPercentage)
                        treshold /= Alpha3; // increase rapidly
                    else if (acceptanceRatio > LargePercentage)
                        treshold *= Alpha2; // decrease slowly
                }
            } while ( (MaxExchanges == 0 ) ? (++CompletedCycles < MaxCycles) : (NrOfExchanges < MaxExchanges) ); // if MaxExchanges is set properly then ignore MaxCycles stopping criterium
            EndResult = lhdXBest; // Set the EndResult
        }
    }
}
